Björn Reineking, 2023-08-05
Model the relative density of animals (also called range distribution or utilisation distribution) as a function of environmental predictors.
We will use the buffalo data set. Install animove metapackage https://github.com/AniMoveCourse/animove_R_package install.packages(“remotes”) remotes::install_github(“AniMoveCourse/animove_R_package”)
This code needs ctmm version 1.1.1 (or newer) remotes::install_github(“ctmm-initiative/ctmm”) Loading packages
library(animove)
library(ctmm)
library(sf)
library(mgcv)
library(mvtnorm)
library(lubridate)
See Kami’s slides for how we got here…
data(buffalo_utm)
data(buffalo_env)
raster::plot(buffalo_env)
raster::plot(raster(buffalo_env, 1))
points(buffalo_utm)
## Loading required package: move
## Loading required package: geosphere
## Loading required package: sp
## Loading required package: raster
##
## Attaching package: 'raster'
## The following object is masked from 'package:nlme':
##
## getData
## The following objects are masked from 'package:ctmm':
##
## projection, projection<-
##
## Attaching package: 'move'
## The following object is masked from 'package:ctmm':
##
## speed
To speed up analyses, we will only work with one individual, Cilla
cilla <- buffalo_utm[["Cilla"]]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)
The first day of Cilla shows some unrealistic movements. We remove those observations
cilla <- cilla[timestamps(cilla) > min(timestamps(cilla)) + days(1), ]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)
Create telemetry object
cilla_telemetry <- as.telemetry(cilla)
## Minimum sampling interval of 3 minutes in Cilla
Fit ctmm model
cilla_guess <- ctmm.guess(cilla_telemetry, CTMM=ctmm(isotropic = TRUE), interactive = FALSE)
cilla_select <- ctmm.select(cilla_telemetry, cilla_guess)
Fit akde
cilla_akde <- akde(cilla_telemetry, cilla_select)
## Default grid size of 3 minutes chosen for bandwidth(...,fast=TRUE).
plot(cilla_akde)
Create named list of rasters
be <- list("elev" = raster(buffalo_env, "elev"),
"slope" = raster(buffalo_env, "slope"),
"var_NDVI" = raster(buffalo_env, "var_NDVI"))
The integrator = “Riemann” option is still in testing, we use it here because it is much faster
cilla_rsf_riemann <- rsf.fit(cilla_telemetry, cilla_akde, R = be, integrator = "Riemann")
## Warning in .local(x, ...): This function is only useful for Raster* objects
## with a longitude/latitude coordinates
## Maximizing likelihood @ n=166035
## Calculating Hessian
summary(cilla_rsf_riemann)
## $name
## [1] "OUF"
##
## $DOF
## mean area diffusion speed
## 8.323012 16.577960 72.507827 71.911902
##
## $CI
## low est high
## var_NDVI (1/var_NDVI) -188.04790701 -72.68967041 42.66856618
## slope (1/slope) -0.37986728 -0.03499134 0.30988461
## elev (1/elev) -0.04693585 -0.01277054 0.02139478
## area (square kilometers) 230.24434776 398.33132029 611.73252692
## τ[position] (days) 4.56389649 7.82595173 13.41956827
## τ[velocity] (minutes) 40.67981657 43.18035547 45.83459945
## speed (kilometers/day) 11.89900607 13.45303449 15.00482017
## diffusion (square kilometers/day) 4.14577993 5.29369812 6.57977920
A suitability map
suitability_riemann <- ctmm::suitability(cilla_rsf_riemann, R=be, grid=crop(be[[1]], extent(cilla) * 2))
raster::plot(suitability_riemann)
Range distribution (includes the ranging behaviour)
my_grid <- crop(be[[1]], extent(cilla) * 2)
agde_cilla <- agde(CTMM = cilla_rsf_riemann, R = be, grid = my_grid)
plot(agde_cilla)
# Plot raster of range distribution
agde_raster <- ctmm::raster(agde_cilla, DF = "PMF")
plot(agde_raster)
Functions to generate quadrature points and predict with the model
rsf_points <- function(x, UD, R = NULL, n = 1e5, k = 1e6, type = "Riemann",
rmax =6*sqrt(UD@CTMM$sigma[1,1]),
interpolation = FALSE) {
# Samples background points from a 2D normal distribution fitted to the relocation data, and extracts environmental
# information from a raster object "R"
# x: telemetry object
# UD: UD object
# R: raster* object
# n: number of background points to sample
# k: weight of presence points
# rmax: maximum distance for Riemann-type integration
# interpolation: do interpolation when sampling the grid
# When type 0´= "MonteCarlo", importance sampling is done
stopifnot(UD@CTMM$isotropic)
stopifnot(type %in% c("Riemann", "MonteCarlo"))
if (type == "Riemann") {
quadrature_pts <- getValues(R)
xy <- as.data.frame(xyFromCell(R, 1:ncell(R)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(R))
xy <- st_transform(xy, CRS(UD@info$projection))
xy <- st_coordinates(xy)
r <- sqrt(((xy[,1] - UD@CTMM$mu[1]))^2 + ((xy[,2] - UD@CTMM$mu[2]))^2)
bg <- data.frame(case_ = 0,
x_ = xy[r<rmax,1], y_ = xy[r<rmax,2], w_ = prod(res(R)), k_ = k)
bg <- cbind(bg, quadrature_pts[r<rmax,])
bg <- sf::st_as_sf(bg, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx <- data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k)
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
xx <- rbind(bg, xx)
} else {
quadrature_pts <- MASS::mvrnorm(n, mu = UD@CTMM$mu, Sigma = UD@CTMM$sigma)
xx <- data.frame(case_ = 0, x_ = quadrature_pts[, 1], y_ = quadrature_pts[, 2], w_ = UD@CTMM$sigma[1,1]/n, k_ = k)
xx <- rbind(xx, data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k
))
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
}
xy <- st_coordinates(xx)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
xx <- cbind(xy, xx)
xx
}
predict_rsf <- function(model, UD, object, include_avail = TRUE, data_crs = UD@info$projection) {
if(include_avail) {
xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
xy <- st_transform(xy, data_crs)
xy <- st_coordinates(xy)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
} else {
xy <- cbind("x_" = rep(0, ncell(object)), "y_" = rep(0, ncell(object)))
}
newdata <- as.data.frame(cbind(xy, getValues(object)))
lambda <- as.numeric(predict(model, newdata, type = "link"))
r <- raster(object, 1)
r[] <- exp(lambda - max(lambda, na.rm = TRUE))
r <- r / sum(getValues(r), na.rm = TRUE)
r
}
Generate quadrature points (“background points”)
set.seed(2)
rsf_cilla_df <- rsf_points(cilla_telemetry, cilla_akde, buffalo_env, interpolation = TRUE)
rsf_cilla_df <- rsf_cilla_df[!is.na(rsf_cilla_df$slope),]
Fit a downweighted Poisson regression
m_rsf_cilla <- glm(case_*k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + slope + var_NDVI,
family = poisson(), data= rsf_cilla_df, weights = w_)
## Warning: glm.fit: fitted rates numerically 0 occurred
Summary of model and confidence intervals for parameter estimates
summary(m_rsf_cilla)
##
## Call:
## glm(formula = case_ * k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev +
## slope + var_NDVI, family = poisson(), data = rsf_cilla_df,
## weights = w_)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.32151 3.56117 -2.898 0.00375 **
## x_ 0.31218 0.36103 0.865 0.38721
## y_ 0.48703 0.44872 1.085 0.27775
## I(-(x_^2 + y_^2)/2) 1.14409 0.27988 4.088 4.35e-05 ***
## elev -0.01275 0.01742 -0.732 0.46426
## slope -0.03507 0.17592 -0.199 0.84201
## var_NDVI -72.71439 58.86105 -1.235 0.21670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1125.8 on 289060 degrees of freedom
## Residual deviance: 1057.5 on 289054 degrees of freedom
## AIC: 1071.5
##
## Number of Fisher Scoring iterations: 24
confint(m_rsf_cilla)
## 2.5 % 97.5 %
## (Intercept) -16.6341489 -2.9214401
## x_ -0.3812531 1.0607780
## y_ -0.3346315 1.4178533
## I(-(x_^2 + y_^2)/2) 0.6772977 1.7773232
## elev -0.0487675 0.0168988
## slope -0.4896424 0.2023228
## var_NDVI -185.7287720 45.8497752
Map of suitability, including the home ranging behaviour
suitability_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2))
raster::plot(suitability_glm)
Map of suitability, without the home ranging behaviour
suitability_no_avail_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(suitability_no_avail_glm)
Comparison of estimates from the two approaches (ctmm and “classic” via glm)
vars <- c("elev", "slope", "var_NDVI")
# Estimates are not identical but similar
cilla_rsf_riemann$beta[vars]
## elev slope var_NDVI
## -0.01277054 -0.03499134 -72.68967041
coef(m_rsf_cilla)[vars]
## elev slope var_NDVI
## -0.01274543 -0.03506537 -72.71439217
# Confidence intervals are of similar width, but location is different
ci_glm <- confint(m_rsf_cilla)
ci_glm[vars, ]
## 2.5 % 97.5 %
## elev -0.0487675 0.0168988
## slope -0.4896424 0.2023228
## var_NDVI -185.7287720 45.8497752
summary(cilla_rsf_riemann)$CI[sapply(vars, function(x) grep(x, rownames(summary(cilla_rsf_riemann)$CI))), ]
## low est high
## elev (1/elev) -0.04693585 -0.01277054 0.02139478
## slope (1/slope) -0.37986728 -0.03499134 0.30988461
## var_NDVI (1/var_NDVI) -188.04790701 -72.68967041 42.66856618
Create maps of suitability based on parameter estimates
M <- getValues(buffalo_env)
eta_glm <- M[,vars] %*% coef(m_rsf_cilla)[vars]
eta_rsf <- M[,vars] %*% cilla_rsf_riemann$beta[vars]
suit_r_glm <- raster(buffalo_env, 1)
suit_r_glm[] <- exp(eta_glm)
suit_r_rsf <- raster(buffalo_env, 1)
suit_r_rsf[] <- exp(eta_rsf)
suit_raster <- raster::stack(suit_r_rsf, suit_r_glm)
names(suit_raster) <- c("ctmm", "classic")
raster::plot(suit_raster)
raster::plot(suit_raster, ext = extent(cilla) * 2)